Problemem, który zostanie poddany analizie w prezentowanym raporcie jest wpływ różnych czynników na rozmiar śledzia oceanicznego wyławianego w Europie. Na przestrzeni ostatnich lat zauważono, że rozmiar ryby stopniowo spada. Zbiór analizowanych danych przedstawia pomiary śledzi oraz warunków, w jakich żyły na przestrzeni ostatnich 60 lat. W ramach połowu jednostki losowo wybierano od 50 do 100 sztuk trzyletnich śledzi.
Podczas przeprowadzania analizy w pierwszym kroku uzupełnione zostały niepełne dane. Następnie w postaci histogramów zaprezentowane zostały rozkłady wartości poszczególnych atrybutów. Przedstawienie atrybutów w postaci macierzy korelacji pozwoliło zobserwować silną zależność między długością śledzi, a temperaturą przy powierzchni wody.
Podczas tworzenia modelu predykcyjnego, który na podstawie zawartych w zbiorze danych cech miał przewidziec rozmiar śledzi wykorzystany algorytm RandomForest. Model został nauczony na treningowym zbiorze danych stanowiącym 80% wszystkich posiadanych danych, a następnie przetestowany na pozostałych 20% danych.
Po dokonaniu analizy możemy zaobserwować, że najbardziej wyraźny wpływ na wielkość wyławianych śledzi ma temperatura przy powierzchni wody (sst) - jej stopniowy wzrost bezpośrednio przyczynił się do spadku średniej wielkości śledzi. Dodatkowo bardzo ważnym czynnikiem okazała się dostępność planktonu - szczególnie widłonogów, której spadek rownież wpłynął na rozmiar śledzi.
library(ggplot2)
library(DT)
library(kableExtra)
library(dplyr)
library(plotly)
library(caret)
library(gridExtra)
library(randomForest)
Kod zapewniający powtarzalność eksperymentów
set.seed(23)
Kod pozwalający na załadowanie danych z pliku csv. Dane zawierają 52583 rekordy opisane 16 atrybutami:
| Kolumna | Opis |
|---|---|
| length | długość złowionego śledzia [cm] |
| cfin1 | dostępność planktonu [zagęszczenie Calanus finmarchicus gat. 1] |
| cfin2 | dostępność planktonu [zagęszczenie Calanus finmarchicus gat. 2] |
| chel1 | dostępność planktonu [zagęszczenie Calanus helgolandicus gat. 1] |
| chel2 | dostępność planktonu [zagęszczenie Calanus helgolandicus gat. 2] |
| lcop1 | dostępność planktonu [zagęszczenie widłonogów gat. 1] |
| lcop2 | dostępność planktonu [zagęszczenie widłonogów gat. 2] |
| fbar | natężenie połowów w regionie [ułamek pozostawionego narybku] |
| recr | roczny narybek [liczba śledzi] |
| cumf | łączne roczne natężenie połowów w regionie [ułamek pozostawionego narybku] |
| totaln | łączna liczba ryb złowionych w ramach połowu [liczba śledzi] |
| sst | temperatura przy powierzchni wody [°C] |
| sal | poziom zasolenia wody [Knudsen ppt] |
| xmonth | miesiąc połowu [numer miesiąca] |
| nao | oscylacja północnoatlantycka [mb] |
raw_data <- read.csv(url("http://www.cs.put.poznan.pl/dbrzezinski/teaching/sphd/sledzie.csv"), na.strings = "?")
kable(head(raw_data, 10))%>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
| X | length | cfin1 | cfin2 | chel1 | chel2 | lcop1 | lcop2 | fbar | recr | cumf | totaln | sst | sal | xmonth | nao |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 23.0 | 0.02778 | 0.27785 | 2.46875 | NA | 2.54787 | 26.35881 | 0.356 | 482831 | 0.3059879 | 267380.8 | 14.30693 | 35.51234 | 7 | 2.8 |
| 1 | 22.5 | 0.02778 | 0.27785 | 2.46875 | 21.43548 | 2.54787 | 26.35881 | 0.356 | 482831 | 0.3059879 | 267380.8 | 14.30693 | 35.51234 | 7 | 2.8 |
| 2 | 25.0 | 0.02778 | 0.27785 | 2.46875 | 21.43548 | 2.54787 | 26.35881 | 0.356 | 482831 | 0.3059879 | 267380.8 | 14.30693 | 35.51234 | 7 | 2.8 |
| 3 | 25.5 | 0.02778 | 0.27785 | 2.46875 | 21.43548 | 2.54787 | 26.35881 | 0.356 | 482831 | 0.3059879 | 267380.8 | 14.30693 | 35.51234 | 7 | 2.8 |
| 4 | 24.0 | 0.02778 | 0.27785 | 2.46875 | 21.43548 | 2.54787 | 26.35881 | 0.356 | 482831 | 0.3059879 | 267380.8 | 14.30693 | 35.51234 | 7 | 2.8 |
| 5 | 22.0 | 0.02778 | 0.27785 | 2.46875 | 21.43548 | 2.54787 | NA | 0.356 | 482831 | 0.3059879 | 267380.8 | 14.30693 | 35.51234 | 7 | 2.8 |
| 6 | 24.0 | 0.02778 | 0.27785 | 2.46875 | 21.43548 | 2.54787 | 26.35881 | 0.356 | 482831 | 0.3059879 | 267380.8 | 14.30693 | 35.51234 | 7 | 2.8 |
| 7 | 23.5 | 0.02778 | 0.27785 | 2.46875 | 21.43548 | 2.54787 | 26.35881 | 0.356 | 482831 | 0.3059879 | 267380.8 | 14.30693 | 35.51234 | 7 | 2.8 |
| 8 | 22.5 | 0.02778 | 0.27785 | 2.46875 | 21.43548 | 2.54787 | 26.35881 | 0.356 | 482831 | 0.3059879 | 267380.8 | 14.30693 | 35.51234 | 7 | 2.8 |
| 9 | 22.5 | 0.02778 | 0.27785 | 2.46875 | 21.43548 | 2.54787 | 26.35881 | 0.356 | 482831 | 0.3059879 | 267380.8 | 14.30693 | 35.51234 | 7 | 2.8 |
Sekcja prezentująca podstawowe statystyki dla zbioru danych:
| Atrybut | Wartość |
|---|---|
| Liczba kolumn | 16 |
| Liczba wierszy | 52582 |
summary(raw_data[,2:16])
## length cfin1 cfin2 chel1
## Min. :19.0 Min. : 0.0000 Min. : 0.0000 Min. : 0.000
## 1st Qu.:24.0 1st Qu.: 0.0000 1st Qu.: 0.2778 1st Qu.: 2.469
## Median :25.5 Median : 0.1111 Median : 0.7012 Median : 5.750
## Mean :25.3 Mean : 0.4458 Mean : 2.0248 Mean :10.006
## 3rd Qu.:26.5 3rd Qu.: 0.3333 3rd Qu.: 1.7936 3rd Qu.:11.500
## Max. :32.5 Max. :37.6667 Max. :19.3958 Max. :75.000
## NA's :1581 NA's :1536 NA's :1555
## chel2 lcop1 lcop2 fbar
## Min. : 5.238 Min. : 0.3074 Min. : 7.849 Min. :0.0680
## 1st Qu.:13.427 1st Qu.: 2.5479 1st Qu.:17.808 1st Qu.:0.2270
## Median :21.673 Median : 7.0000 Median :24.859 Median :0.3320
## Mean :21.221 Mean : 12.8108 Mean :28.419 Mean :0.3304
## 3rd Qu.:27.193 3rd Qu.: 21.2315 3rd Qu.:37.232 3rd Qu.:0.4560
## Max. :57.706 Max. :115.5833 Max. :68.736 Max. :0.8490
## NA's :1556 NA's :1653 NA's :1591
## recr cumf totaln sst
## Min. : 140515 Min. :0.06833 Min. : 144137 Min. :12.77
## 1st Qu.: 360061 1st Qu.:0.14809 1st Qu.: 306068 1st Qu.:13.60
## Median : 421391 Median :0.23191 Median : 539558 Median :13.86
## Mean : 520367 Mean :0.22981 Mean : 514973 Mean :13.87
## 3rd Qu.: 724151 3rd Qu.:0.29803 3rd Qu.: 730351 3rd Qu.:14.16
## Max. :1565890 Max. :0.39801 Max. :1015595 Max. :14.73
## NA's :1584
## sal xmonth nao
## Min. :35.40 Min. : 1.000 Min. :-4.89000
## 1st Qu.:35.51 1st Qu.: 5.000 1st Qu.:-1.89000
## Median :35.51 Median : 8.000 Median : 0.20000
## Mean :35.51 Mean : 7.258 Mean :-0.09236
## 3rd Qu.:35.52 3rd Qu.: 9.000 3rd Qu.: 1.63000
## Max. :35.61 Max. :12.000 Max. : 5.08000
##
no_na_data <- sum(complete.cases(raw_data))
na_data <- nrow((raw_data)) - no_na_data
Możemy zauważyć, że w badanym zbiorze danych aż 10094 rekordów posiada brakujące dane. Niestety usunięcie takiej ilości rekordów z badanego zbioru mogłoby wpłynąć negatywnie na jakość analizy dlatego konieczne jest uzupełnienie brakujących wartości. Analizując przedstawione statystyki możemy zaobserwować, że jedynie niektóre kolumny posiadają brakujące wartości. Są to dane dotyczące dostępności planktonu oraz temperatury przy powierzchni wody. To właśnie dane w tych cechach należy uzupełnić.
Przyglądając się danym możemy zauważyć, że braki występują w sąsiedztwie danych kompletnych. Wynika to prawdopodobnie z faktu, że sąsiadujące rekordy pochodzą z tego samego połowu, a braki w danych wynikają z nieuwagi podczas ich zbierania.
Skrypt uzupełniający dane pobierany jest z pliku “ZED_grouping.R”.
Dla każdego rekordu z brakującą wartością cechy pobierane jest najbliższe 6 kompletnych rekordów. Na podstawie pobranej grupy wyliczana jest mediana brakującej wartości, która następnie wstawiana jest jako wartość badanej cechy w niekompletnym rekordzie. Czynność powtarzana jest dla każdej cechy zawierającej brakujące dane.
# Załadowanie skryptu uzupełniającego dane
if(!exists("fillMissingValues", mode="function")) source("ZED_grouping.R")
filled_data <- fillMissingValues(raw_data)
filled_data
ggplot(filled_data, mapping = aes(x = cfin1)) + geom_histogram(binwidth = 0.01)
ggplot(filled_data, mapping = aes(x = cfin2)) + geom_histogram(binwidth = 0.01)
ggplot(filled_data, mapping = aes(x = chel1)) + geom_histogram(binwidth = 0.01)
ggplot(filled_data, mapping = aes(x = chel2)) + geom_histogram(binwidth = 0.01)
ggplot(filled_data, mapping = aes(x = lcop1)) + geom_histogram(binwidth = 0.01)
ggplot(filled_data, mapping = aes(x = lcop2)) + geom_histogram(binwidth = 0.01)
ggplot(filled_data, mapping = aes(x = fbar)) + geom_histogram(binwidth = 0.01)
ggplot(filled_data, mapping = aes(x = recr)) + geom_histogram(binwidth = 60)
ggplot(filled_data, mapping = aes(x = cumf)) + geom_histogram(binwidth = 0.01)
ggplot(filled_data, mapping = aes(x = totaln)) + geom_histogram(binwidth = 60)
ggplot(filled_data, mapping = aes(x = sst)) + geom_histogram(binwidth = 0.01)
ggplot(filled_data, mapping = aes(x = sal)) + geom_histogram(binwidth = 0.005)
ggplot(filled_data, mapping = aes(x = xmonth)) + geom_histogram(binwidth = 0.01)
ggplot(filled_data, mapping = aes(x = nao)) + geom_histogram(binwidth = 0.05)
Korelacja między atrybutami została przedstawiona w postaci macierzy korelacji.
kable(cor(filled_data), digits = 2)%>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"))%>%
row_spec(0,angle = -45)
| X | length | cfin1 | cfin2 | chel1 | chel2 | lcop1 | lcop2 | fbar | recr | cumf | totaln | sst | sal | xmonth | nao | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| X | 1.00 | -0.34 | -0.15 | 0.06 | -0.16 | 0.05 | -0.23 | 0.04 | 0.09 | 0.00 | 0.23 | -0.36 | 0.36 | -0.06 | 0.00 | 0.41 |
| length | -0.34 | 1.00 | 0.08 | 0.10 | 0.22 | -0.01 | 0.24 | 0.05 | 0.25 | -0.01 | 0.01 | 0.10 | -0.45 | 0.03 | 0.01 | -0.26 |
| cfin1 | -0.15 | 0.08 | 1.00 | 0.15 | 0.09 | 0.20 | 0.12 | 0.21 | -0.06 | 0.12 | -0.05 | 0.13 | 0.01 | 0.13 | 0.01 | 0.01 |
| cfin2 | 0.06 | 0.10 | 0.15 | 1.00 | 0.00 | 0.31 | -0.04 | 0.65 | 0.15 | -0.10 | 0.34 | -0.22 | -0.24 | -0.08 | 0.01 | -0.01 |
| chel1 | -0.16 | 0.22 | 0.09 | 0.00 | 1.00 | 0.29 | 0.95 | 0.25 | 0.16 | -0.05 | 0.07 | 0.17 | -0.22 | -0.15 | 0.05 | -0.51 |
| chel2 | 0.05 | -0.01 | 0.20 | 0.31 | 0.29 | 1.00 | 0.17 | 0.88 | 0.03 | 0.00 | 0.26 | -0.38 | 0.01 | -0.22 | 0.07 | -0.06 |
| lcop1 | -0.23 | 0.24 | 0.12 | -0.04 | 0.95 | 0.17 | 1.00 | 0.15 | 0.10 | 0.00 | -0.01 | 0.27 | -0.26 | -0.10 | 0.03 | -0.55 |
| lcop2 | 0.04 | 0.05 | 0.21 | 0.65 | 0.25 | 0.88 | 0.15 | 1.00 | 0.05 | 0.00 | 0.29 | -0.30 | -0.12 | -0.18 | 0.06 | -0.04 |
| fbar | 0.09 | 0.25 | -0.06 | 0.15 | 0.16 | 0.03 | 0.10 | 0.05 | 1.00 | -0.24 | 0.82 | -0.51 | -0.18 | 0.04 | 0.01 | 0.07 |
| recr | 0.00 | -0.01 | 0.12 | -0.10 | -0.05 | 0.00 | 0.00 | 0.00 | -0.24 | 1.00 | -0.26 | 0.37 | -0.20 | 0.28 | 0.02 | 0.09 |
| cumf | 0.23 | 0.01 | -0.05 | 0.34 | 0.07 | 0.26 | -0.01 | 0.29 | 0.82 | -0.26 | 1.00 | -0.71 | 0.03 | -0.10 | 0.03 | 0.23 |
| totaln | -0.36 | 0.10 | 0.13 | -0.22 | 0.17 | -0.38 | 0.27 | -0.30 | -0.51 | 0.37 | -0.71 | 1.00 | -0.29 | 0.15 | -0.03 | -0.39 |
| sst | 0.36 | -0.45 | 0.01 | -0.24 | -0.22 | 0.01 | -0.26 | -0.12 | -0.18 | -0.20 | 0.03 | -0.29 | 1.00 | 0.01 | -0.01 | 0.51 |
| sal | -0.06 | 0.03 | 0.13 | -0.08 | -0.15 | -0.22 | -0.10 | -0.18 | 0.04 | 0.28 | -0.10 | 0.15 | 0.01 | 1.00 | -0.02 | 0.13 |
| xmonth | 0.00 | 0.01 | 0.01 | 0.01 | 0.05 | 0.07 | 0.03 | 0.06 | 0.01 | 0.02 | 0.03 | -0.03 | -0.01 | -0.02 | 1.00 | 0.00 |
| nao | 0.41 | -0.26 | 0.01 | -0.01 | -0.51 | -0.06 | -0.55 | -0.04 | 0.07 | 0.09 | 0.23 | -0.39 | 0.51 | 0.13 | 0.00 | 1.00 |
Kod umożliwiający analizę rozmiaru śledzi w czasie.
W celu analizy rozmiaru rozmiaru śledzi w czasie, na podstawie informacji o rocznym narybku oraz wiedzy, że dane występują w porządku chronologicznym utworzona została dodatkowa kolumna “year”.
if(!exists("add_year", mode="function")) source("ZED_grouping.R")
mean_sample <- filled_data[,c("recr","length")]
recr_groups <- aggregate(mean_sample, by = list(mean_sample$recr), FUN = mean)
year_data <- add_year(filled_data, recr_groups)
aggreagated_years <- aggregate(year_data, by= list(year_data$year), FUN = mean)%>%select(year, length, lcop2, lcop1, chel2, sst)
interactive_plot <- ggplot(aggreagated_years, mapping = aes(x = year, y = length))+geom_smooth()
ggplotly(interactive_plot)
W następnym kroku podjęta została próba utworzenie regresora, który przewidzi rozmiar śledzi. Dane zostały podzielone na zbiór uczący oraz testowy w proporcji 80%/20%. Wykorzystany został algorytm RandomForest optymalizujący model według miary RMSE.
sig_attr_selection <- filled_data%>%select(length, cfin1, cfin2, chel1, chel2, lcop1, lcop2, fbar, recr, cumf, totaln, sst, sal, nao)
inTraining <- createDataPartition(y = sig_attr_selection$length, p=.8, list=FALSE)
training_data <- sig_attr_selection[inTraining,]
test_data <- sig_attr_selection[-inTraining,]
rfGrid <- expand.grid(mtry = 1:10)
ctrl <- trainControl(
method = "repeatedcv",
number = 5,
repeats = 5)
fit <- train(length ~ .,
data = training_data,
method = 'rf',
trControl = ctrl,
metric = "RMSE",
tuneGrid=rfGrid,
importance = TRUE,
ntree=10)
ggplot(fit) + theme_bw()
rfClasses <-predict(fit, newdata = test_data[-1])
rmse <- RMSE(pred = rfClasses, obs = test_data$length)
r2 <- R2(pred = rfClasses, obs = test_data$length)
Algorytm osiągnął najlepszy rezultat dla parametru “mtry” równego 3.
| Miara | Wartość |
|---|---|
| RMSE | 1.1825109 |
| R2 | 0.4916135 |
Na podstawie wyuczonego modelu oszacowane zostały wagi każdego z atrybutów zbioru danych.
attr_importance <- varImp(fit, scale = FALSE)
ggplot(attr_importance)
Możemy przyjrzeć się zmiane najbardziej istotnych cech w czasie.
lcop2_plot <- ggplot(filled_data, mapping = aes(x = X, y = lcop2))+geom_smooth()+ggtitle("Dostępność planktonu [zagęszczenie widłonogów gat. 2]")
ggplotly(lcop2_plot)
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
chel2_plot <- ggplot(filled_data, mapping = aes(x = X, y = chel2))+geom_smooth()+ggtitle("Dostępność planktonu [zagęszczenie Calanus helgolandicus gat. 2]")
ggplotly(chel2_plot)
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
lcop1_plot <- ggplot(filled_data, mapping = aes(x = X, y = lcop1))+geom_smooth()+ggtitle("Dostępność planktonu [zagęszczenie widłonogów gat. 1]")
ggplotly(lcop1_plot)
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
fbar_plot <- ggplot(filled_data, mapping = aes(x = X, y = fbar))+geom_smooth()+ggtitle("Ułamek pozostawionego narybku")
ggplotly(fbar_plot)
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
sst_plot <- ggplot(filled_data, mapping = aes(x = X, y = sst))+geom_smooth()+ggtitle("Temperatura przy powierzchni wody [°C]")
ggplotly(sst_plot)
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
Oczywiście wpływ na wielkość wyławianych śledzi miało wiele czynników, jednak wśród najbardziej wpływowych wyróżnić możemy dostępność planktonu oraz temperaturę przy powierzchni wody.